TO APPEAR IN THE ASTROPHYSICAL JOURNAL 

Preprint typeset using LATf^X style emulateapj v. 11/12/01 



MODIFICATION OF PROJECTED VELOCITY POWER SPECTRA BY DENSITY 
INHOMOGENEITIES IN COMPRESSIBLE SUPERSONIC TURBULENCE 

Christopher M. Brunt 1 ' 2 ' 3 & Mordecai-Mark Mac Low 4 

To appear in the Astrophysical Journal 

ABSTRACT 

The scaling of velocity fluctuation Sv as a function of spatial scale I in molecular clouds can be measured 
from size-linewidth relations, principal component analysis, or line centroid variation. Differing values 
of the power law index of the scaling relation {Sv 2 } 1 / 2 oc P 3D in three dimensions (3D) are given by 

£f) \ these different methods: the first two give 7315 ~ 0.5, while line centroid analysis gives 73,0 ~ 0. This 

discrepancy has previously not been fully appreciated, as the variation of projected velocity line centroid 
fluctuations ((Sv 2 ^ 1 / 2 oc P 2D ) is indeed described, in 2D, by j 2 d ~ 0.5. However, if projection smoothing 

04 , is accounted for, this implies, in 3D, that 73c w 0. We suggest that a resolution of this discrepancy 

■ can be achieved by accounting for the effect of density inhomogeneity on the observed 72D obtained 
from velocity line centroid analysis. Numerical simulations of compressible turbulence are used to show 
that the effect of density inhomogeneity statistically but not identically reverses the effect of projection 
smoothing in the case of driven turbulence so that velocity line centroid analysis does indeed predict 

0^ I that 72D ~ 73D ~ 0.5. For decaying turbulence, the effect of density inhomogeneity on the velocity line 

centroids diminishes with time such that at late times, 72D ~ 73D + 0.5 due to projection smoothing 
alone. Deprojecting the observed line centroid statistics thus requires some knowledge of the state of the 
flow. This information can be inferred from the spectral slope of the column density power spectrum and 
I ■ a measure of the standard deviation in column density relative to the mean. Using our numerical results 

\^ I we can restore consistency between line centroid analysis, principal component analysis and size-linewidth 

■^j- , relations, and we derive 73c «s 0.5, corresponding to shock-dominated (Burgers) turbulence, which also 

t— I ■ describes the simulations of driven turbulence at scales where numerical dissipation is negligible. We 

t-H I find that this consistency requires that molecular clouds are continually driven on large scales or are 

CO . only recently formed. 

Subject headings: ISM: clouds — ISM: kinematics and dynamics — ISM: molecules — methods: 

■ statistical — radio lines: ISM — turbulence 
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we estimate the observed power spectrum of line centroid 
velocity fluctuations in the NGC 7129 molecular cloud, 
and interpret the measured power spectrum exponent and 
previously measured line centroid structure function ex- 
ponents, accounting both for projection smoothing and 
density inhomogeneity. 

2. VELOCITY LINE CENTROID ANALYSIS : THEORY 

There have been a number of correlation studies on pro- 
jected line centroid velocities in molecular clouds (Kleiner 
& Dickman 1985; Kitamura et al. 1987; Hobson 1992; Mi- 
esch & Bally 1994; Ossenkopf & Mac Low 2002). Most 
of the studies have proceeded in direct space, relying in 
large part on the formalism established by Kleiner & Dick- 
man (1985) and Dickman & Kleiner (1985). Here we use 
the related Fourier space representation, and measure line 
centroid velocity fluctuations via their power spectrum. 
The power spectrum (P(k) as a function of wavenum- 
ber k) of velocity fluctuations is: P(k) oc k~ KND , where 
knd = ^Ind + N in TV dimensions 5 . The above relation- 
ship between knd and -f^o applies to the case of exact 
power law power spectra. In general, the power spectrum 
and structure function of fluctuations are integrally to each 
other (the structure function is a simple transformation of 
the autocorrelation function which in turn is the inverse 
Fourier transform of the power spectrum) and defining the 
statistical relationship between them with single indices is 
difficult. Our observational data does have well-defined 
power law behavior over the observed range of scales, so 
the above simple prescription is valid. 

We use the power spectrum for three reasons: (1) 
in the numerical fields, increasing separations ultimately 
decrease due to field periodicity resulting in artificially 
smaller jnd for a given Knd (see Brunt et al. 2003); 
(2) the numerical fields have power spectra that are not 
power laws over all wavenumbers due to the scales at which 
the turbulence was driven and the effects of numerical 
dissipation- we use a Fourier space scheme that avoids 
problems arising from this; (3) for the observational data, 
corrections for instrumental noise and beam-smearing are 
easily achieved in Fourier space-any direct space beam- 
smearing corrections must in general proceed via Fourier 
space anyway. We thus characterize all measurements in 
terms of power spectra, but note the direct space equiva- 
lent according to the simple kj^o~Jnd given above. 

In the power spectrum representation, the Kolmogorov 
case has k^d = 11/3 and the Burgers case has k^d = 4. 
Previously, the power spectrum approach has been used 
by Kitamura et al. (1987), and, more recently, by Miville- 
Deschenes et al. (2003(b)) for H I 21 cm line emission. 

The true projected velocity line centroid field is : 

vic ue (x,y) = v z{x,y,zi) / N pix , (1) 

i=l 

where v z (x,y,z) is the line of sight velocity component 
and N pix is the number of spatial pixels along each line 
of sight (this is replaced by a similar integral form for the 
continuous case; our numerical and observational data are 
pixelized, so we will use the pixelized equations here). 

5 Note that for kjvd < N, 7jvd saturates to ~ zero, while for kmd > 
et al. (2003) for a demonstration of the latter case. 



Consider first an unweighted projection of the velocity 
field over the line of sight (z) axis. Generation of a velocity 
centroid map, by projection over the z direction, simply 
selects the wavevectors for which k z = 0, since all others 
are oscillatory along the z direction and will average out 
(disregarding inhomogencous sampling). For an isotropic 
power law power spectrum in 3D, the projected velocity 
line centroid 2D power spectrum has the same spectral 
slope as it does in 3D; i.e. «2D = k-sd (see e.g. Stutzki 
et al. 1998; Miville-Deschenes et al. 2003; Brunt et al. 
2003); more precisely, the observed (projected) 2D power 
spectrum is identical to the k z = plane in the 3D power 
spectrum. Only the transverse component of the velocity 
field is projected into 2D (e.g. Kitamura et al. 1987); in 
other words : the longitudinal contribution to v z has no 
amplitude on the k z = plane. 

For the case of a power law power spectrum, while the 
spectral slope is unchanged on going from 3D to projected 
2D, the direct space statistical properties of the field are 
modified. The variation of root mean square velocity fluc- 
tuation ((Sv 2 ) 1 ' 2 ) with spatial scale (£) is : 

{8V 2 ) 1 ' 2 oc ^n D -n)/2 K pND ^ ^ 

where N is the dimension on which the field is defined 
(here, N = 2 or 3) and where v, for the 2D case, is under- 
stood to be the projected line centroid field, v\ r c ue . The 
criterion k 3 d — k 2 d leads to -f2D = 73D + 0.5; this is 
projection smoothing. For the Kolmogorov (1941) case 
(k,3D = K-2D = 11/3 ; 73£> = 1/3) the projected index is 
J2D = 5/6 (von Hoerner 1951; see O'Dcll & Castancda 
1987). For Burgers (1974) compressible turbulence, we 
have k 3 d — — 4, 73D = 0.5, and -f2D = 1 due to 
projection smoothing. 

An observed 72D = 0.5 thus implies 7313 = 0. Only in 
the case that the line of sight depth over which the obser- 
vational tracer exists is much less than the total transverse 
scale of emission observed on the sky (Miville-Deschenes 
et al. 2003; see also O'Dell & Castaneda 1987) is it true 
that 72_d ~ 73r>; i-e. this requires thin, sheet-like clouds, 
oriented face-on towards us. 

The density weighted projected velocity line centroid 
field is : 

v 7c 9t ( x >y) = ^2n(x,y, Zi ) v z (x,y,z l ) / ^ n ( x iV> z i)i 
i=i i=i 

(3) 

where n(x, y, z) is the density field. For the case of uni- 
formly excited, optically thin gas, Equation 3 is the appro- 
priate form for predicting the velocity line centroid field. 
As recognized by Kleiner & Dickman (1985), the form of 
Equation 3 suggests that the projected "velocity line cen- 
troid" field is more akin to a measure of the line-of-sight 
momentum of the medium; however, the normalization by 
the projected column density complicates this assignment 
(dimensionally, the line centroid field has units of velocity). 
For uniform density (e.g. the case of fully incompressible 
turbulence) the density field amplitude is normalized out 
of the line centroid field. More generally, for inhomoge- 
neous density, the mean projected density field is normal- 
ized out, but variations in density along the line of sight 

N + 2, 7jvd saturates to unity - see Stutzki et al. (1998) and see Brunt 
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may induce additional structure in the measured velocity 
line centroid field. 

Previously established results of relevance here are : (1) 
Miville-Deschenes et al. (2003) showed that density in- 
homogeneity does not affect velocity line centroids in the 
case of uncorrelated density and velocity fields (thus pre- 
serving = 73_d + 0.5). However Miville-Deschenes et 
al. (2003) used stochastic density and velocity fields (in- 
cluding some degree of non-physical correlation between 
density and velocity). (2) Brunt et al. (2003) found that 
density inhomogeneity does induce additional small-scale 
structure on the velocity line centroid field in the case 
of physically correlated density and velocity (resulting in 
72D < 73D + 0.5). They also found that the correlations 
between density and velocity are important (for a given 
amplitude of density fluctuations). Brunt et al. (2003) 
used "spectrally-modified" numerically-simulated density 
and velocity fields (see Lazarian et al. 2001) in which the 
density- velocity coupling is also modified to some extent. 
(3) Ossenkopf & Mac Low (2002) examined numerical sim- 
ulations of turbulence, and measured 73 £> in 3D and 72 d 
from the projected line centroid velocity field of the sim- 
ulations (including density- weighting) . They found that 
73D ~ 12D ~ 0.5; they also noted that projection smooth- 
ing does not predict this behavior. This is an indication 
that the effect of density weighting counters the effects of 
projection smoothing. (4) Lazarian & Esquivel (2003, sub- 
mitted) have inferred that density inhomogeneity docs not 
affect the power spectrum of velocity line centroid fluctu- 
ations. 

Observationally, the density weighted line centroid ve- 
locity field may be approximated by use of an optically 
thin tracer : 

N chan N chan 

v°c s {x,y) = T ( x ,y, v zi) v Zi (x,y) / ^ T(x,y,v Zi ) 

i=l i=l 

(4), 

where T(x,y,v Zi ) is the observed radiation temperature, 
v Zi (x, y) is the velocity of the i th channel at position x, y, 
and N c han is the number of channels along each line of 
sight. 

3. NUMERICAL DATA 

We use simulations of randomly driven hydrodynami- 
cal (HD) turbulence and magnetohydrodynamical (MHD) 
turbulence (Mac Low 1999), and decaying hydrodynami- 
cal turbulence (Mac Low et al. 1998). These simulations 
were performed with the astrophysical MHD code ZEUS- 
3D 6 (Clarke 1994), a 3D version of the code described 
by Stone & Norman (1992a, b) using second-order advec- 
tion (van Leer 1977), that evolves magnetic fields using 
constrained transport (Evans & Hawley 1988), modified 
by upwinding along shear Alfven characteristics (Hawley 
& Stone 1995). The code uses a von Neumann artificial 
viscosity to spread shocks out to thicknesses of three or 
four zones in order to prevent numerical instability, but 
contains no other explicit dissipation or resistivity. Struc- 
tures with sizes close to the grid resolution are subject to 
the usual numerical dissipation, however. 

The simulations were performed on a 3D, uniform, 
Cartesian grid with periodic boundary conditions in ev- 

6 Available from the Laboratory for Computational Astrophysics at http: 



ery direction. All length scales are in fractions of the cube 
size; i.e. all wavevectors run from k = 1 to k = N p i X /2 for 
a simulation of N p i X pixels in each dimension. An isother- 
mal equation of state was used in the computations, with 
sound speed chosen to be c s =0.1 in arbitrary units. The 
initial density and, where applicable, magnetic field were 
both initialized uniformly on the grid, with the initial den- 
sity po = 1 and the initial field parallel to the z-axis. 

The turbulent flow is initialized with velocity pertur- 
bations drawn from a Gaussian random field determined 
by its power distribution in Fourier space, following the 
usual procedure. As discussed in detail in Mac Low et al. 
(1998), the decaying turbulence runs were initialized with 
a flat spectrum with power from kd = 1 to kd = 8 be- 
cause that will decay quickly to a turbulent state. Note 
that the dimcnsionless wavenumber kd counts the number 
of driving wavelengths in the box. A fixed pattern of 
Gaussian fluctuations drawn from a field with power only 
in a narrow band of wavenumbers around some value kd 
offers a very simple approximation to driving by mecha- 
nisms that act on that scale. For driven turbulence, such 
a fixed pattern was normalized to produce a set of per- 
turbations 5v(x,y, z). At every time step a velocity field 
Sv(x,y,z) — A5v was added to the velocity v, with the 
amplitude A now chosen to maintain constant kinetic en- 
ergy input rate, as described by Mac Low (1999). 

A summary of the models is given in Table 1. 

4. OBSERVATIONAL DATA 

Our observational data are 12 CO (J=l-0) and 
13 CO (J=l-0) spectral line imaging observations acquired 
with the Five College Radio Astronomy Observatory 
(FCRAO) 14 m telescope during the Spring season, 2003. 
These data form part of a large (~ 200 square degree) 
molecular line survey (in progress) within the ongoing 
multiwavelength Canadian Galactic Plane Survey (CGPS; 
Taylor et al. 2003). 

The FCRAO 14 m telescope was equipped with the SE- 
QUOIA focal plane array (Erickson et al. 1999), which 
was recently expanded to 32 pixels. The array is rectan- 
gular (4x4 pixels) with two sets of 16 pixels in differ- 
ent polarizations having the same pointing centers. Due 
to the expanded physical size of the array in the dewar, 
it is no longer possible for the dewar to rotate and so 
maintain a constant orientation with the celestial coor- 
dinate system (our mapping grid here is in Galactic l,b 
coordinates). To allow mapping in celestial coordinates, 
FCRAO has developed an on-the-fly (OTF) mapping sys- 
tem in which the telescope is driven across the map area 
and spectra are dumped on the desired grid. One pixel 
near the center of the array (central pixel hereafter) is se- 
lected to define the origin of the map coordinates. The 
other pixels traverse the map along paths determined by 
the pixel separation in the array and the current orien- 
tation of the celestial coordinate system to the telescopic 
(azimuth-elevation; az-el) system. If the central pixel is 
scanned so that it Nyquist samples the map area, then the 
combination of data from all the other pixels results in a 
highly (but irregularly) oversampled map. In order to map 
large areas, the sampling rate of the central pixel can be 
relaxed slightly. Our data are obtained by dumping the 

: / /cosmos. ucsd.edu/lca/ 
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correlators every 22.5" (1/2 beamwidth at the 12 CO fre- 
quency) along the scan direction during a total scan length 
(in I) of 0.5°. The b coordinate of the central pixel is then 
incremented by 33.75" and the next scan in I is taken. The 
pixel separation in the array is 88.56". In the case that 
the I, b system is aligned with the az-el system, then after 
two scans, the first row of the array is located at 66.45" in 
b above its original position, or 22.11" below the original 
b coordinate of the second row. The next scan places the 
first row of the array 11.64" above the first scan of the sec- 
ond row of the array and 12.11" below the second scan of 
the second row. For subsequent scans, a similar progres- 
sion of sampling ensues. Apart from at the extreme map 
edges, a sampling finer than the Nyquist rate is achieved, 
but the sampling is slightly irregular. More generally, the 
I, b system is rotated with respect to the az-el system and 
an even finer sampling is achieved. An OFF (reference) 
measurement is taken between every two scans of the ar- 
ray (i.e. OFF, scan out and back, OFF, scan out and 
back, etc), and a calibration measurement is taken after 
every 4 "out and back" sequences (« every 9-10 minutes). 
The data are calibrated to the T\ scale using the standard 
chopper wheel method (Kutner & Ulich 1981). Following 
recording of the data, the spectra are converted onto a reg- 
ular grid of 22.5" pixel spacing using the FCRAO otftool 
software. 

The 12 CO and 13 CO spectral lines were acquired simul- 
taneously using the recently added dual channel correla- 
tors, which provide a total bandwidth of 50 MHz, with 
1024 channels at each frequency. The central velocity of 
the spectra is vi sr = -45 km s _1 , resulting in a total veloc- 
ity coverage of V[ sr « -110 km s _1 to vi sr « +20 km s _1 
with a channel spacing of ~ 49 kHz (« 0.13 km s _1 ). The 
spatial resolution of the data (beam FWHM) is 45" for 
12 CO and 46" for 13 CO. 

The otftool software was configured to remove first order 
baselines derived from the spectra in the velocity intervals: 
-110 to -95 km s _1 and +7 to +20 km s" 1 . Spectra from 
individual pixels in the array contributing to a single po- 
sition in the output grid were weighting by the reciprocal 
of their noise variance derived from the above velocity in- 
tervals. Following initial gridding, we removed additional 
(« sinusoidal) baseline structure arising from a standing 
wave between the subreflector and dewar using Fourier 
methods. First, we defined the "spectral window" between 
-20 km s -1 and +5 km s _1 which contains the detectable 
emission in the region of interest here. Then we interpo- 
lated over this window using a first order baseline derived 
from the adjacent 100 channels on cither side of the spec- 
tral window. The emission (+ noise) in the spectral win- 
dow at each spatial position was replaced by the interpo- 
lated line at that position and the resulting spectrum was 
converted to Fourier space using a Fast Fourier Transform 
(FFT). The first two wavevectors, of wavelength 50 MHz 
and 25 MHz (« 130 and 65 km s _1 respectively) were 
extracted from the FFT'ed spectrum, converted back to 
direct space and then subtracted from the original spec- 
trum. 

The resulting sensitivities (ltr, TJ) of the processed data 
are ~ 0.5 K and - 0.25 K per channel for 12 CO and 13 CO 
respectively. 



5. VELOCITY LINE CENTROID ANALYSIS : NUMERICAL 
RESULTS 

5.1. Overview 

The numerical fields may be used to evaluate the role of 
density weighting on the power spectrum of the line cen- 
troid velocity field. We calculated v\ r c ue via Equation 1 
and v™ c 9t via Equation 3 for all simulations. For the MHD 
models we examined line centroid fields projected over two 
axes: parallel and perpendicular to the mean magnetic 
field direction. As the simulated velocity fields have pe- 
riodic boundary conditions, the k z = plane in the 3D 
power spectrum of v z is identical to the 2D power spec- 
trum for u'£ ue . 

First, wc examine the driven simulations. Example ve- 
locity line centroid fields are shown in Figure 1. It is im- 
mediately evident that the inclusion of density-weighting 
results in the insertion of additional small-scale structure 
into the line centroid fields. 

To quantify the effect seen in Figure 1, we calculated 
the power spectra (P(k) as function of wavenumber k) of 
all the line centroid fields, using a FFT. Figure 2(a) shows 
representative results for HC8 model (kd — 7-8) . There is 
a turnover in the power spectra at or near kd- For k < kd, 
the true and weighted line centroid power spectra are quite 
similar, while for k > kd the power spectrum of v^" c 9t is 
shallower. We can quantify the effect of density weighting 
quite precisely by examining the difference (in log space) 
between the power spectra of v™ gt and v^ ue (in real space 
this corresponds to measuring the spectral index of the ra- 
tio P v w 9 t /P v tr U e). We characterize the ratio of power spec- 
tra by a scaling exponent 5k, where : P v ^gt /P v t ruc cx k~ SK 
(i.e log(P v n, g t)-log(P v t^i) = const + 5k log(kj). This pro- 
cedure avoids potential problems arising from non-power- 
law projected velocity spectra. 

The ratio P v ^gt /P v trv, c for the HC8 model is shown in 

I c I c 

Figure 2(b). There is little modification to the true line 
centroid power spectrum for spatial scales above the driv- 
ing scale (k = 7-8); F ' w 9 t / 1 P v ± r ^ does not vary much with 

k at small k. This is because there are no major density 
enhancements on these large scales (Ballesteros-Paredes 

6 Mac Low 2002). For k > k d , P v w g t / 'P v tme is a well de- 
veloped power law, with some deviation at very high k. 
We fitted the exponent 5k = -1.17 ± 0.06 over the range 

7 < k < 32 (i.e. from the largest driving scale to a wave- 
length of 4 pixels). The largest scale at which dissipation 
effects are first seen in these simulations is ~10 zones or 
k w 12.8 (Ossenkopf & Mac Low 2002) but we do not see 
any spectral breaks in the power spectrum ratio before 
k = 32; we will examine the effects of this in greater detail 
in Section 5.2. 

The other driven models gave similar results. Figure 3 
shows P v ^gtjP v tTuc for all driven models. These spectra 
were fitted between the largest driving scale and k = 32 
for all models to derive 5k, which in all cases was negative 
and ~ -1 with variations of up to ~ ± 0.4. The measured 
5k are listed in Table 1. There is no obvious systematic 
difference between the HD and MHD cases, nor any obvi- 
ous variation with rms Mach number, M, for the driven 
runs. There is some evidence that 5k is more negative 
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when the line of sight is parallel to the mean magnetic 
field direction, but this does not universally apply. Tak- 
ing the observed range of Sk as a natural variation, we 
derive a mean (unweighted) Sk = 0.96 ± 0.22 from all 
the driven models. In other words, inclusion of density 
weighting results in a shallower spectral slope, with index 
~1 less than the true line centroid power spectrum index. 
In direct space, such a modification would thus correspond 
to 72D ~ 73D rather than ^ 2 d = 73 d + 0.5. (In general, 
72D = 73D + 0.5 +5k/2.) 

If the driven turbulence conditions are representative of 
the conditions within molecular clouds, the scaling expo- 
nents measured from line centroid structure functions are 
approximately equal to the true, 3D exponents. We stress 
that the "cancellation" of projection smoothing by small 
scale line centroid structure induced by density inhomo- 
geneity is statistical in nature, and the effect of density 
weighting does not of course identically reverse the effects 
of projection smoothing. 

The same procedure to derive Sk was then applied to 
the simulations of decaying turbulence. Figure 4 shows 
the plots of P wgt /P v tT US evolving as time increases. We 

ic ' c 

have fitted Sk over the range 10 < k < 64 for all fields 
(motivated again by the lack of clear breaks in the power 
spectrum ratio) and these are listed in Table 1. This shows 
that as the turbulence decays, the effect of density inho- 
mogcncity on the line centroid field is reduced, due to the 
lessened density contrast. If molecular clouds are in such 
a state, then projection smoothing only applies to the ob- 
served line centroid field, and 7 2 _d ~ 73/5 + 0.5. Over the 
course of one sound crossing time, the density contrast 
is reduced sufficiently, even for initially highly supersonic 
turbulence, that Sk becomes very small, due to weaker 
shocks in decaying turbulence relative to those in driven 
turbulence (Smith, Mac Low, & Zuev 2000; Smith, Mac 
Low, & Heitsch 2000); in other words, if the turbulence is 
continually driven, or has only recently entered a decaying 
state, then we expect Sk « -1, but this evolves quickly to 
Sk — > over a sound crossing time. Figure 5 summarizes 
the evolution of Sk after the driving is turned off. 

5.2. Regime of Validity 

As noted above, we are fitting Sk over a range of k that 
includes the dissipation region (k greater than <~ 12.8), 
motivated by the lack of distinct breaks in the power spec- 
trum ratio. To test whether dissipativc effects could induce 
systematic biases in our results, we repeated the fits while 
restricting the range to a maximum k = 12. This was 
done only for the driven models for which kd — 1-2 and 
kd = 3-4 in order that sufficient data was still available 
for the fits. The resulting measurements of Sk are plot- 
ted against the measurements of Sk fitted to k — 32 in 
Figure 6. For this sample, we find that the mean Sk are 
-0.91 ± 0.17 (fitted to k = 32) and -0.87 ± 0.19 (fitted to 
k = 12). These are not significantly different, so our visual 
conclusion that there are no breaks in the power spectrum 
ratios is good. 

The advantage of using a relative measure, Sk, in the 
above analysis is that it circumvents potential problems 
arising from small curvatures in the intrinsic velocity 
power spectra. We now examine actual measurements of 
the spectral indices of the projected velocity power spec- 



tra. These absolute indices are needed to examine the 
range of validity of Sk as derived above. Due to the cur- 
vature in the power spectra, likely arising from dissipative 
effects at higher k, the absolute spectral slopes depend 
on the range in k over which they are derived. In what 
follows, we label the spectral slope of the unweighted pro- 
jected velocity field as K2D,u and the spectral slope of the 
density- weighted projected velocity field as k 2 d,w These 
are related by k 2 d, w -k 2 d, u = Sk. 

Example line centroid power spectra for the HC2, HC4 
models, driven at kd = 1-2 and 3-4 respectively, are shown 
in Figure 7. These are fitted over two ranges in k: up to 
k = 32 (upper panels) and up to k = 12 (lower panels); 
in both cases the smallest k used to make the fit is taken 
as the smaller of the driving wavenumbers. When the fit 
is made to k = 32, the spectra are "steep" (k2d,u ~ 4.5) 
but when the fit is made to k = 12, the indices are in 
good agreement with the value expected for Burgers tur- 
bulence (k2d,u ~ 4). Numerical dissipation will cause just 
the sort of steepening of the spectrum (less structure at 
small scale) that is seen in Figure 7. 

We have derived k 2 d,u, k 2 d,w and Sk for the driven 
runs over the two ranges in k given above. Figure 8(a) 
shows the relation between k 2 d,u and k 2 d, w for all driven 
runs fitted to A: = 32; the corresponding variation of Sk 
with k 2 dm is shown in Figure 8(b). Figures 8(c) and 8(d) 
show the same relations obtained when the fits are made to 
k = 12 for the subset of runs with kd = 1—2 and kd = 3-4. 
In all these Figures, the solid line denotes Sk = -0.96 de- 
rived in the preceding Section, and the dashed lines mark 
uncertainties of ±0.22 around this value. 

These tests show that our estimated mean Sk = 
0.96 ± 0.22 is not strongly affected by the fitting range 
in k but there is some uncertainty in the absolute range 
of k 2 d,u and k 2 d^ w to which this value applies. Conserva- 
tively, the fits to k = 12 should be used, as these are free 
of dissipative effects. Using only the fields for which such 
fits can be made, we find that Sk = -0.87 ± 0.19, and this 
applies to k 2 dm ~ 3.5-4.2 and k 2 d, w ~ 2.7-3.4. 

However, we also note that the lack of dependence of Sk 
on the fitting range in k is indirect evidence that a single 
value of Sk is applicable at all k 2 e> m , but this interpreta- 
tion should be examined in more detail in future. In what 
follows we shall assume that Sk = -0.96 ± 0.22 holds for 
all k 2 d, u and k 2 d, w - 

5.3. Origin of the Spectral Modification 
5.3.1. Effect of Density Spectral Index 

Lazarian & Esquivel (2003) have suggested that (in our 
terminology) Sk may have a larger amplitude for "shallow" 
density field power spectra. The density field power spec- 
trum slope may be estimated from the power spectrum 
of the projected column density field, P c d(k), using the 
same arguments as presented in § 2 (see also Lazarian & 
Pogosyan 2000; Stutzki et al. 1998). If Sk is related to the 
spectral slope, k C( j, of the column density power spectrum 
(where P c d(k) cx k~ Kcd ) then this would indeed provide a 
means of inferring Sk observationally. In other words, the 
observed column density power spectrum slope is used as 
a probe of the (relative) density fluctuation spectrum as a 
function of scale within the medium. 

To investigate this, we produced column density maps 
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for all the simulations (including projection along two axes 
- parallel and perpendicular to the magnetic field - for the 
MHD fields). We then fitted K c d for each field, using the 
same range in wavenumber that was used to derive 5k. 
The fitted values of K c d are given in Table 1. (The col- 
umn density power spectra are also slightly curved, so k c< i 
would be lower if fitted to k = 12; thus the K c d axis is 
uncertain by a ~linear shift as in Figure 8 for k 2 d,u-) 

In Figure 9 we plot the variation of 5k with K c d- There 
is a clear variation of 5k with K c d, demonstrated by the 
decaying turbulence simulations. As time proceeds, K c d 
increases (i.e. preferential reduction in small scale density 
contrast) and the amplitude of 5k diminishes (i.e. less- 
ened effect of density inhomogeneity on the velocity line 
centroid field). For the driven simulations, the K c d are 
typically ~ 3. These are in reasonable accord with obser- 
vational estimates of K c d ~ 2.8 for molecular clouds; e.g. 
Bcnsch, Stutzki, & Ossenkopf 2001. For our K c d ~ 3, 5k 
varies around the mean value of -1. Thus, observational 
estimates of K c d suggest that accounting for density in- 
homogeneity is very important for velocity line centroid 
analysis of molecular clouds, implying (with 5k « -1) that 
72D ~ 73D counter to the expectations based on projec- 
tion smoothing alone. This provides an explanation for 
the 72D ~ 73D result observed in these simulations by Os- 
senkopf & Mac Low (2002). Lazarian & Pogosyan do not 
provide a physical motivation for whether the medium is 
density dominated (i.e. "shallow density spectrum" hav- 
ing K c d < 3) or not: it is simply a question of what K c d 
is, regardless of why. In our scenario, this regime is the 
driven turbulence regime. 

5.3.2. Amplitude of Density Fluctuations 

A note of caution regarding the use of K c d to infer 5k 
should be made: K c d is a non-unique measure of the den- 
sity field. To investigate this more deeply, we analyzed 
modified versions of the HE2 field. First, we randomized 
the Fourier space phases of the density and velocity fields 
to produce new fields with the same power spectra as the 
original fields. This modification has two effects : (1) it 
removes the density-velocity correlations, (2) it drastically 
reduces the amplitude of density fluctuations. Following 
phase-randomization of the density field, we found that 
we had to recenter the mean density due to the creation 
of negative densities by the phase-randomization process. 
After renormalization back to a mean density of unity, 
the phase-randomized density field had a standard devi- 
ation of 0.21 and range of lxl0~ 6 to 2.02. For compari- 
son, the original density field, normalized to mean unity, 
had a standard deviation of 2.53 and a range of 1.2xl0 -6 
to 92.1. The logarithm of the phase-randomized density 
field had a standard deviation of 0.22, compared to the 
standard deviation of the logarithm of the original density 
field of 2.49. Using the phase-randomized density fields 
and velocity fields we derived 5k — -0.18 ± 0.03. This 
is substantially lower than the original 5k = -1.04 ± 0.05 
obtained with the original HE2 density and velocity fields. 
This shows that the "shallowness" of the density power 
spectrum, as measured by K c d, does not uniquely allow 
inference of 5k. 

We conducted a second test, this time randomizing the 
phases of the logarithm of density, and then exponentiating 



the phase-randomized log(density) field. This procedure 
again removes the density-velocity correlations but does 
not lead to the creation of negative densities, and allows 
a greater amplitude of density contrast. We found that 
the resulting K c d ~ 2.41 for the density field produced 
by this method was slightly shallower than the original 
K c d ~ 2.87 for the HE2 model. The standard deviation 
of log (density) was the same for the phase-randomized 
density field and the original density field (2.49). Using 
the exponentiated phase-randomized log(density) field in 
conjunction with the phase-randomized velocity field, we 
found 5k = -1.16 ± 0.05. 

These two tests together suggest that it is the amplitude 
of density fluctuations, measured as the standard deviation 
of log(density) that is controlling the value of 5k. To make 
sure that this is the case, and that density- velocity correla- 
tion does not play a major role, we took the original veloc- 
ity field and the original density field, shifted by 64 pixels 
from its original position. For a shift along the linc-of- 
sight, we derived 5k = -0.90 ± 0.04 and for a shift trans- 
verse to the line-of-sight we derived 5k = -0.97 ± 0.06. 
These are slightly lower in magnitude than the original 
5k = -1.04 ± 0.05, but the difference is marginally signif- 
icant. 

We thus conclude that the amplitude of density contrast 
is most important in determining 5k. To examine this di- 
rectly, we measured the standard deviation of log (density), 
(J ;n( / 9)i f° r all the numerical fields (see Table 1) and we 
plotted this versus 5k (Figure 10(a)) and against K c d Fig- 
ure 10(b)). This indicates that there is a fairly good cor- 
relation between 5k and o~i n ( p \, and that with the assump- 
tion that the numerical fields are a better representation 
of real molecular clouds than the phase-randomized fields, 
that Kcd may be used as a surrogate for cr;„( p ). Wc find 
that the driven simulations are characterized by K c d ~ 3 
and (T;„(p) exceeding ~ unity, but that there is no obvious 
dependence of 5k on K c d or er;„( p ) within this regime. 

Ideally, one would like to measure (T;„( p ) to infer 5k but 
it is of course impossible to measure cr ; „( p ) directly How- 
ever, it may be possible to infer a hl ( p ^ indirectly using 
the column density field. Figure 11(a) compares the vari- 
ation in the standard deviation of log(column density), 
o'iniN)! to G in(p)- The dashed line is fitted : cr^^) = - 
0.02 ± 0.08 + (0.33 ± 0.06)ct Mp) as a guide. The values 
of cri„(7v) are listed in Table 1. Unfortunately, it is also 
impossible to use o^at) observationally, as the small N 
are dominated by noise. Instead, we consider simply the 
standard deviation of density, c p / Po , and column density 
°~n/n i where both fields are normalized their mean val- 
ues, po and N 0} respectively. The dependence of cr N / No 
on c p / Po is shown in Figure 11(b) along with the fitted 
line, a N/No = -0.06 ± 0.07 + (0.34 ± 0.04)cr p/po . These 
results show that the variation in column density may be 
used to infer the variation in (3D) density, but the dy- 
namic range in N is compressed due to line-of-sight aver- 
aging (<Tjv/jv ~ °V/po/3)- However, in general, it should 
be recognized that, at any finite resolution, the mean den- 
sity and column density are well-estimated, but their stan- 
dard deviations are continuously variable with resolution 
in the respective dimensions. It may be possible to at- 
tempt to construct a resolution-independent measure to 
compare 2D and 3D standard deviation but this is beyond 
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the scope of the current work. 

In Figure 12 we plot combinations of o- p / po , (Jn/n , 5k, 
and K cc i to gauge the utility of different measures. In sum- 
mary, some guidelines concerning 5k may be established 
observationally by measuring k c( i and <Jn/n q - Note that a 
measurement of K c d and <7jv/jv may also be used to dis- 
tinguish the phase-randomized fields from the numerically 
simulated fields : the phase-randomized density field has 
a much lower <j n / Nq = 0.065 compared with that of the 
numerical field {<tn/n = 0.786) for the same K c d = 2.87. 
Some caution comparing <j n / Nq with Figure 12 is needed 
if the observations are at significantly different spatial dy- 
namic range than the simulations used here. 

5.4. Comparison with Previous Results 

Our results appear inconsistent with two previous stud- 
ies. Brunt et al. (2003) found that 5k w -0.4 from the 
simulation used in their paper (Mach number w 1). The 
density field used in Brunt et al. (2003) had <7;„( p ) = 0.45. 
Comparison with Figure 10(a) shows that these numbers 
are not discrepant with the current results. The conclu- 
sion of Brunt et al. (2003) that, for a given amplitude of 
density fluctuations, the correlations between density and 
velocity are important was based on the a single simula- 
tion and is outweighed by our new results here. Secondly, 
our results are inconsistent with the results (5k « 0) found 
by Lazarian & Esquivel (2003) for their simulation (Mach 
number w 2.5). This may be due to the solenoidal driving 
in their simulations (Cho & Lazarian 2002); certainly this 
discrepancy deserves further attention. 

5.5. Resolution Tests 

Finally, we examined the effect of limited spatial dy- 
namic range on our results. The driven numerical fields 
implicitly contain realizations at different effective resolu- 
tions (i.e. due to the varying driving scales). These show 
no obvious trend of 5k with driving wavelength; however 
there is a large amount of scatter in the derived 5k that 
may hide such a trend. In addition we have explored the 
effect of limiting the fits to 5k to different ranges of k. We 
conducted a further, more direct, tests using realizations 
of model D (decaying turbulence) at sizes of 128 3 and 64 3 
and a realization of model HC2 (driven turbulence) at 256 3 
resolution. These realizations utilize the same parameter- 
ization but are not identical in detail (i.e. they were not 
initialized with the same random number seeds and so are 
subject to the natural scatter in 5k evidenced in Table 1). 

Snapshots of model D at 128 3 and 64 3 sizes are available 
at timesteps of 0.1, 0.15 and 0.3 sound crossing times. We 
conducted the same analysis to derive 5k from these low- 
ered resolution simulations. At 128 3 size, we find 5k = - 
0.78 ± 0.08, -0.57 ± 0.04 , and -0.22 ± 0.03 for times 
0.1, 0.15 and 0.3 respectively. At 64 3 size, we find 5k = - 
0.71 ± 0.21, -0.33 ± 0.22 , and -0.11 ± 0.25 for times 0.1, 
0.15 and 0.3 respectively As the resolution is increased, 
the 5k are converging to the values of 5k =-0.76 ± 0.03- 
0.59 ± 0.02, and -0.32 ± 0.02 found at 256 3 size at the 
same timesteps. The driven model HC2, with higher Mach 
number, showed a somewhat more pronounced shift, from 
5k = -0.72 ± 0.05 at 128 3 to 5k = -1.01 ± 0.03 at 256 3 . In 
sum, this resolution study suggests that the magnitude of 
5k is decreased slightly by limited spatial dynamic range. 



Therefore, our derived 5k magnitudes should probably be 
considered to be lower limits (i.e. the true, high resolution, 
values may be a little more negative than found here). 

6. VELOCITY LINE CENTROID ANALYSIS : 
OBSERVATIONAL RESULTS 

6.1. Overview 

An overview of the CGPS CO data in the vicinity of 
NGC 7129 is given in Figure 13. We have selected a re- 
gion (64 x 64 pixels, 32 x 32 beams, 4096 spectra) for 
velocity line centroid analysis, indicated by the superim- 
posed square in Figure 13. This region is located away 
from the main star forming area (it contains no IRAS point 
sources), and has simple, singly-peaked spectral lines. To 
derive velocity line centroid maps, we utilized Gaussian 
fitting to the spectral lines rather than a direct projection 
via Equation 4, as Gaussian fitting was found to be less 
noisy (see below). Example fits to both the 12 CO lines 
and 13 CO lines are shown in Figure 14. These spectra 
are taken from a 16 x 16 pixel block at the center of the 
analysis region and are block-averaged over 4x4 pixels. 
The 12 CO lines have roughly twice the peak temperature 
of the 13 CO lines, resulting in comparable signal-to- noise 
for each tracer. 

Fits to the 13 CO lines could not be made for the en- 
tirety of the analysis region, due to weak or absent emis- 
sion. Use of the 13 CO line centroid field for analysis would 
then require values for the missing data to be assumed be- 
fore calculation of the power spectrum. A possible choice 
would be to set the missing data to the mean value of the 
measurable data. However, a better procedure is to create 
a "hybrid" line centroid field, using the 12 CO centroid ve- 
locities as the best estimate of the missing 13 CO centroid 
velocities. We thus masked off all pixels in the analysis 
region for which the peak temperature of the 13 CO line 
was less than 1.4 K Tjj[ (i.e. ~ 5<r), and replaced the ab- 
sent or noisy 13 CO measurements with the 12 CO centroid 
velocities. For the hybrid field, 40% of the line centroid 
velocities are from 13 CO and 60% are from 12 CO. 

For the lines of sight where we have good fits to both 
13 CO and 12 CO, we compared the derived line centroid 
velocities, as shown in Figure 15(a). The two measure- 
ments are very similar over the entire variation of centroid 
velocity (~ -11 to -9.5 km s _1 ) but there is a slight av- 
erage offset between the two, of -0.075 km s _1 (~ 1/2 
channel width) as seen in Figure 15(b). This can arise 
from greater saturation on one side of the 12 CO line pro- 
files, or from depth-dependent abundance variations in the 
presence of a line-of-sight velocity gradient. The effects of 
opacity on our measured 12 CO velocity line centroid map 
are thus very small. With the offset understood, the width 
of the distribution in Figure 15(b) (~ ± 0.1 km s _1 ) can 
be otherwise attributed to noise on the measurements. To 
measure the noise on the line centroid fields, we took the 
fitted Gaussian line profiles as the best estimate of the true 
line profiles and added actual instrumental noise from the 
same spatial positions, taken from a signal-free area of the 
baseline region. We then carried out the same Gaussian 
fitting procedure to derive the noise distribution. For com- 
parison, we derived the direct line centroid field from the 
"re-noised" fields using Equation 4. The estimated errors 
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on the 12 CO line centroid field are ± 0.093 km s _1 (Gaus- 
sian fitting) and ± 0.122 km s _1 (using Equation 4). The 
Gaussian-fitting procedure is less noisy, and as the lines 
in the original data are well-described by Gaussians (Fig- 
ure 14), we used the line centroids derived by Gaussian fit- 
ting for the power spectrum analysis below. The resulting 
velocity line centroid fields derived from Gaussian fitting 
are shown in Figure 16. The observed fields are qualita- 
tively more similar to the simulations driven at large scale 
(the kd — 1-2 case in Figure 1) than to those driven at 
small scale. 

6.2. Large Scale Gradients and Edge Effects 

Prior to the full analysis, we investigated the effect of 
filtering out the "large-scale" line centroid velocity fluc- 
tuations. Whether this should be done or not is a mat- 
ter of some controversy (Miesch & Bally 1994; Osscnkopf 
& Mac Low 2002). We estimated the large-scale gradi- 
ents by fitting a quadratic surface to the 12 CO line cen- 
troid map. Then we subtracted this from the original line 
centroid map, and calculated the power spectrum before 
and after subtraction, using FFTs. We found that only 
the low frequency (k = 1) power is strongly affected by 
this procedure, as may have been expected. The mea- 
sured spectral slopes, between k — 2 and k = 13 (see 
below) are k 2 o = 3.38 ± 0.10 (before subtraction) and 
k 2 d = 3.42 ± 0.09 (after subtraction). Thus the spectrum 
steepens slightly (but insignificantly) after subtraction of 
the large-scale gradients, most likely due to the suppres- 
sion of high frequency structure resulting from mismatch 
in the line centroids at the field edges (i.e. the FFT rou- 
tine assumes periodic fields). As there is little difference 
between the line centroid power spectrum before and after 
subtraction of the large-scale gradients, and it is not clear 
whether they should be subtracted, we do not employ the 
subtraction in the analysis below. 

We also examined the effect of apodization (Miville- 
Deschenes, Lagache & Puget 2002). We centered the 12 CO 
line centroid field to zero mean, then applied a cosine ta- 
per to a strip 2 pixels wide around the edge of the field. 
The edges of the line centroid field now match at zero 
velocity. For the apodized field, we applied the same ta- 
per function to the noise estimates, to ensure that the 
noise floor is still well estimated. At large k both the 
apodized line centroid power spectrum and the apodized 
noise floor are reduced, resulting in comparable signal-to- 
noise before and after apodization. After subtraction of 
the respective noise floors (see below) the fitted spectral 
slopes are k 2 d — 3.70 ± 0.13 (before apodization) and 
= 3.76 ± 0.12 (after apodization). There is little 
to distinguish these measurements over the measurement 
uncertainty, so we do not utilize apodization in the anal- 
ysis below. As the apodization routine simply imposes 
arbitrary smooth structure in the line centroid field, it is 
not surprising that the measured power spectrum slope 
steepens slightly (i.e. reduction of small scale power). Al- 
ternatively, the untapered field contains unphysical sharp 
velocity structure at the field edges, so the actual noise- 
floor subtracted spectral slope most likely lies between the 
two slope estimates above. 

6.3. Correction for Noise and Beam Smearing 



The raw measured indices stated above (k 2 d ~ 3.4) 
are not the actual true measured indices, as we have not 
yet accounted for the effects of measurement noise and 
beam smearing on the measured power spectra. The to- 
tal power spectrum in the presence of noise is, simply, 
Ptrue{k) + Pnoise(k) (where P true is the power spectrum 
that would be observed in the absence of noise) assuming 
independence of the signal and noise. The noise power 
spectrum must thus be subtracted from the raw mea- 
sured power spectrum prior to analysis. The finite size 
of the telescope beam induces correlations between closely 
spaced pixels; for a beam pattern B{x, y) in direct space 
(i.e. roughly a Gaussian form with 45" FWHM), the ob- 
served line centroid map (neglecting noise for now) is : 
v ic S — v°^ s ' true <8> B where (g> denotes convolution and 
v obs,tme j g ^ e K ne centroid field that would be observed 
without beam smearing. The beam smeared power spec- 
trum is then P abs(k) = P [ absque (fc)Pe(A:) where Pb(/c) is 

lc V lc 

the power spectrum of the beam pattern. The final output 
power spectrum, corrected for noise and beam smearing is: 

Pout P ut(k) = (P <s (k) - P nolse (k)) / P B (k). (5) 

Note that the division out of the beam pattern is applied 
after the noise floor subtraction as the noises in our data 
are spatially independent (in contrast to the measured sig- 
nal). 

The noise power spectra can be derived from the "re- 
noised" fields discussed above; we did this for both the 
12 CO line centroid field and the hybrid line centroid field. 
We assumed that a Gaussian beam pattern, B(x,y) of 
45" FWHM applied to both the 12 CO line centroid field 
and the hybrid line centroid field. A summary of the rele- 
vant terms in Equation 5 is provided in Figure 17. Output 
power spectra were then derived via Equation 5 and these 
are shown in Figure 18. The resulting scaling range covers 
k = 2 to k = 13 limited at low k by a possible turnover 
and the uncertainty of large-scale gradient subtraction, 
and limited at high k by the amplification of measure- 
ment noise by beam pattern division. The resulting fits 
to the range k = 2 to k — 13 are k 2 d = 3.11 ± 0.10 and 
k 2D = 2.93 ± 0.13 for the 12 CO line centroid field and the 
hybrid line centroid field respectively. These are jointly 
consistent with k 2 d ~ 3, giving 72D ~ 0.5. 

The range of the fit appears deceptively large in Fig- 
ure 18; the spatial dynamic range achieved is only 6.5, 
between wavelengths of 32 and 4.9 spatial pixels (^2.5 
beamwidths). Similar spatial dynamic range was achieved 
by Miesch & Bally (1994), who measured a typical value of 
72D ~ 0.43. Around three orders of magnitude in spatial 
dynamic range was achieved for the Polaris Flare molecu- 
lar cloud by Ossenkopf & Mac Low (2002) by combining 
three data sets at different angular resolutions and spatial 
scales; they also measured 72D ~ 0.5. 

6.4. Interpretation 

If we interpret the measured index of k 2 d ~ 3 
(72U ~ 0.5) allowing for projection smoothing but without 
taking into account the effect of density weighting then 
the inferred 3D spectral index is k^d w 3. We would thus 
obtain 73^ m (via Equation 2) or : no dependence of 
"turbulent" velocity dispersion on scale. If correction for 
projection smoothing is applied to the indices measured by 
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Miesch & Bally (1994) and Ossenkopf & Mac Low (2002), 
this again gives 73c w 0. A similar result (i.e. 730 w 0) 
was found by O'Dell & Castaneda (1987) for H II region 
turbulence, after correcting for projection smoothing. The 
correction for projection smoothing, however, has been the 
exception rather than the rule for velocity line centroid 
analysis and most authors have taken their projected 72D 
as being equivalent to the intrinsic 730- This is not cor- 
rect, as the selection of the k z — plane of the 3D line 
of sight velocity power spectrum by projection onto the 
2D plane of the sky is purely a geometrical effect : pro- 
jection smoothing thus always applies (unless the cloud 
is actually two dimensional - see below), being modified 
only by the effects of density inhomogeneity (and radia- 
tive transfer). Hobson (1992) interprets (e.g.) k 2 d = 3 
(i.e. 72D = 0.5) as being applicable to shock-dominated 
turbulence (73c = 0.5). Miesch & Bally (1994) directly 
compare their 72D to theoretical predictions of 730 with 
no modification. Lazarian & Esquivcl (2003) derive from 
their numerical data (in our terminology) 5k = 0, but then 
interpret the 72D ~ 0.43 of Miesch & Bally (1994) as mean- 
ing also that 73 d « 0.43. 

Our numerical results here show that the (statistical) 
equivalence of 72D and 73^ is actually valid in the case 
of driven turbulence but only due to the accidental can- 
cellation of two competing factors : projection smoothing, 
which suppresses small scale velocity line centroid struc- 
ture, and density inhomogeneity, which creates small scale 
velocity line centroid structure. 

The interpretation of existing results from line centroid 
analysis applied to molecular clouds in general will de- 
pend on the inferred, or assumed, state of turbulence. If 
molecular clouds are in an evolved state of decaying tur- 
bulence, then projection smoothing only applies, and we 
would infer 73D « for all clouds studied to date. If, 
however, the clouds are still being strongly driven, or have 
only recently formed, then the 5k s=s -1 induced by density 
inhomogeneity should apply, and we thus can find a bet- 
ter consistency between line centroid analysis and other 
measures of 73c- To show this, we applied principal com- 
ponent analysis (PCA) to the analysis region according to 
the prescription given in Brunt & Heycr (2002) to derive a 
sequence of coupled characteristic scales, 5v and £, in ve- 
locity and space respectively. The results of this, yielding 
slightly better than an order of magnitude in spatial dy- 
namic range, are shown in Figure 19. The 5v and I pairs 
were fitted as 5v oc l a . We derive a = 0.58 ± 0.01 ( 13 CO) 
and a = 0.59 ± 0.02 ( 12 CO), resulting in j 3D =0.44 ± 0.01 
( 13 CO) and 7 3D =0.46 ± 0.04 ( 12 CO) according to the 
Brunt (1999) calibration. Some caution is necessary here, 
though, on two counts. First, a is more directly depen- 
dent on the first order velocity fluctuation spectrum, not 
the second order spectrum that is traced by the power 
spectrum of the line centroid field; second, estimates of a 
are uncertain (around the mean calibration line) by ~0.1 
(Brunt et al. 2003). Neither of these factors can reconcile 
the PCA results with 730 « however. If, as suggested by 
PCA in the limit of weak intermittency (as also supported 
by the lack of strong exponential tails in the line profiles), 
73D ~ 0.45 then we infer that 5k ~ -0.9 would be needed 
to reconcile the line centroid results with PCA. This value 
of 5k is consistent with 5k — -0.96 ± 0.22 derived from 



the numerical simulations of driven turbulence. 

To determine the observationally-motivated value of 5k, 
without recourse to PCA, we examined the power spectra 
of the integrated intensity fields for both isotopes (it would 
be very difficult to attempt this with a hybrid integrated 
intensity field). Here, we assume that the integrated inten- 
sity field provides a good estimate of the column density 
field. To do this, we estimated the noise floor by integrat- 
ing signal-free regions of the baseline over the same num- 
ber of channels as was used to generate integrated intensity 
maps. The measured "column density" power spectra were 
then corrected for the noise floor and beam pattern (Equa- 
tion 5) and the resulting output power spectra are shown 
in Figure 20. We have fitted K c d over the range k = 2 
to k = 13 (the same as used for the velocity line centroid 
power spectra), to estimate K c d « 2.5 consistently between 
the two isotopes. This is in general agreement with typical 
K c d measured in molecular clouds (Bensch et al. 2001) but 
slightly smaller than can be accounted for by our numer- 
ical fields. Using Figure 9 as a rough guide, we interpret 
the measured K c d as motivating 5k « -1. In doing this, 
we are considering the numerically simulated fields to be a 
more realistic representation of real molecular clouds than 
the phase-randomized versions. We do not see a turnover 
at small k in the integrated intensity power spectra, nor 
in the velocity line centroid power spectra, except perhaps 
at k = 1 which is quite uncertain. Our observed power 
spectra are thus more in accord with the numerical mod- 
els driven at small wavenumber (see also Ossenkopf & Mac 
Low 2002). Finally, we used the 13 CO data to provide an 
estimate of ct/v/at , by assuming that the 13 CO is optically 
thin. To do this, we created an integrated intensity im- 
age of the emission and an integrated intensity image of 
instrumental noise over the same number of channels (64); 
the data were left in T4 units, since the absolute scal- 
ing is unimportant. The mean integrated intensity of the 
data is : I = 1.034 ±0.004 K km s" 1 . The total (sig- 
nal+noise) standard deviation of the integrated intensity 
image is o-j, to t = 0.661. The contribution of the standard 
deviation of the noise is : oi^oise — 0.253 K km s _1 . 
Thus we measure 07 signal/ i a — 0.591 ±0.004 which we 
take as being m o~M/N a - This is clearly inconsistent with 
o~N/N a = 0.065 seen in the phase-randomized density field 
(obtained with K c d = 2.87) and referral to Figure 12(d) 
will show a good consistency between estimates of K c d and 
cjv/JVo f° r the observational data compared to the numeri- 
cal data. The effective spatial dynamic range in the simu- 
lations for computing ctn/No wm ne between the nominal 
spatial dynamic range (128 per axis) and the most strin- 
gent interpretation of the maximum dissipation scale (~ 
12.8 per axis). These compare reasonably well with the 
data, at a spatial dynamic range of 32 beams (64 pixels) 
per axis. We interpret the consistency found between K c d 
and <7n/n as motivating the use of the numerically simu- 
lated fields to infer 5k s=a -1. 

Our results regarding the effect of density inhomogene- 
ity indicate that the measured k 2 d ~ 3 is thus shallower 
than the true k 3 d by approximately 5k ~ -1, and we then 
infer k 3 d ~ 4 (c.f. Figure 8(c)) and thus -f 3 o « 0.5. 
In more detail, we derive (using 5k = -0.96 ± 0.22) 
k 3D = 4.07 ± 0.24 ( 12 CO) and k 3D = 3.89 ± 0.26 (hy- 
brid). These convert to : 731, = 0.54 ± 0.12 ( 12 CO) and 
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73D = 0.45 ± 0.13 (hybrid). These are in good agreement 
with 73£> 0.45 determined from PCA, and are consistent 
with compressible, Burgers turbulence. The uncertainty in 
the inferred 73^1 is large, due mainly to the intrinsic vari- 
ability of 5 k (~ ± 0.22). We thus cannot entirely rule out 
the Kolmogorov case k 3 d ~ 3.67, 73D ~ 0.33, from veloc- 
ity line centroid analysis alone, particularly if we include 
the extra uncertainties arising from subtraction of large 
scale gradients and edge effects (estimated to be ~ ± 0.1 
in k 2 d)- 

We interpret the conditions required for consistency 
between velocity line centroid analysis, principal compo- 
nent analysis and size-linewidth analysis as providing in- 
direct support for molecular clouds that have only re- 
cently formed and/or arc continually driven on large scales 
(Elmegreen 1993; Ballcsteros-Paredes, Vazquez-Semadeni, 
& Scalo 1999). 

Projection smoothing does not apply in the case that the 
line of sight depth over which the observational tracer ex- 
ists is much less than the total transverse scale of emission 
observed on the sky, requiring sheet-like clouds, oriented 
face-on towards us. For a face-on sheet, the effect of den- 
sity weighting would also vanish as there is no longer any 
line of sight variability (see Equation 3 in the limit that 
Npi X — ► 1 as a representation of this regime) . While this 
solution is quite unsatisfactory in general, it represents a 
limiting case : no projection smoothing, no density weight- 
ing effects, resulting in the direct measurement of intrinsic 
statistical properties; i.e. the field is actually two dimen- 
sional. Some modified version of this limiting case may 
be operating in our numerical fields; i.e. in the case that 
a single density enhancement (not necessarily confined to 
a thin sheet, but having some spatial coherence) is domi- 
nant along the line of sight, then this will approximate the 
Npix — > 1 limit of Equation 3. 

In general, for the full 3D (isotropic) driven case, our 
numerical results show that projection smoothing is coun- 
tered by the effects of density inhomogeneity, again result- 
ing in recovery of intrinsic statistical scaling properties, 
but in a less direct way. 

7. SUMMARY 

We have examined numerical simulations of supersonic 
turbulence in molecular clouds in order to understand the 
effects of projection smoothing and density inhomogeneity 
on observable velocity line centroid fields (vi c (x,y)). We 
find that the observed index, -f 2 D, describing the variation 
of projected line centroid velocities ((Svf^ 1 ^ 2 oc p 2D ) is re- 
lated to the intrinsic index, 73^, (where {Sv 2 ) 1 / 2 cx P 3D ) 



via the following relation : 

72D = 73D + 0.5 + Sk/2, 

where the increase of 0.5 is a direct consequence of pro- 
jection smoothing, and 5k is an empirically determined 
correction for the effects of density inhomogeneity. 

For driven turbulence, both with and without magnetic 
fields, 5k rs -1 (with variations ~ ± 0.22), thus implying 
12D ~ 73D due to the accidental cancellation of the effects 
of projection smoothing and small-scale structure induced 
by density inhomogeneity on the line centroid field. For 
decaying turbulence, 5k tends to ~0 at late times, so that 
12D ~ 73D + 0.5 due entirely to projection smoothing. 

The variations in 5k observed in high Mach number 
driven turbulence have no obvious dependence on Mach 
number, M. As M — >0 (the incompressible case) one 
should also expect that 5k — >0, as is indeed seen at late 
times in the decaying runs. We have inferred that 5k is de- 
termined mostly by the amplitude of density fluctuations. 
The scatter in 5k in the supersonic regime emphasizes the 
need for continued high spatial dynamic range surveying 
of the molecular ISM in order to provide a statistically 
useful sample of fields. 

Deprojection of observed scaling relations in velocity 
line centroid fields requires some knowledge of the state 
of the flow. This information can be inferred from the 
spectral slope of the column density power spectrum and 
a measure of the standard deviation in column density 
relative to the mean {<Jn/n )- Our numerical results show 
that 73_d is now consistently estimated at ^0.5 using ve- 
locity line centroid analysis and using principal component 
analysis and size-linewidth relations. This consistency re- 
quires that molecular clouds are continually driven or only 
recently formed. 
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Table 1 



Numerical Models : Parameters and Measurements 



Model 




HD /MHD /Decay 


k 


b 

d 


Af 


va/c/ 


1 t° 




)K 




i 




! 




<p) 


&ln(N) 


HA8 


128 


HD 


7- 


-8 


1.9 





1 


-0.84 


± 


0.05 


4.10 


± 


0.12 





.71 


0.17 


HC2 


128 


HD 


1 


-2 


7.4 





1 


-0.72 


± 


0.05 


3.30 


± 


0.13 


1 


.58 


0.71 


HC4 


128 


HD 


3 


4 


5.3 





1 


-0.98 


± 


0.03 


3.22 


± 


0.09 


1 


.55 


0.48 


HC8 


128 


HD 


7- 


-8 


4.1 





1 


-1.17 


± 


0.06 


3.36 


± 


0.08 


1 


.34 


0.32 


HE2 


128 


HD 


1- 


-2 


15.0 





0.475 


-1.04 


± 


0.05 


2.87 


± 


0.10 


2 


.45 


0.74 


HE4 


128 


HD 


3- 


-4 


12.0 





0.875 


-1.11 


± 


0.06 


3.02 


± 


0.09 


2 


.23 


0.63 


HE8 


128 


HD 


7- 


-8 


8.7 





1 


-1.28 


± 


0.11 


3.14 


± 


0.09 


2 


.07 


0.40 


MC4X v ± 


128 


MHD 


3- 


-4 


5.3 


10 


0.1 


-0.91 


± 


0.12 


3.04 


± 


0.09 


1 


.68 


0.56 


MC4X U|| 


128 


MHD 


3- 


-4 


5.3 


10 


0.1 


-1.21 


± 


0.05 


2.57 


± 


0.09 


1 


.68 


0.46 


MC45 v± 


128 


MHD 


3 


-4 


4.8 


5 


0.2 


-0.66 


± 


0.09 


3.03 


± 


0.13 


1 


.58 


0.45 


MC45 v n 


128 


MHD 


3 


-4 


4.8 


5 


0.2 


-1.06 


± 


0.05 


2.93 


± 


0.11 


1 


.58 


0.48 


MC41 v ± 
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MHD 


3 


-4 


4.7 


1 


0.5 


-0.89 


± 


0.05 


3.18 


± 


0.08 


1 


.28 


0.36 


MC41 u|| 
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MHD 


3 


-4 


4.7 


1 
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± 
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3.20 


± 
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1 


.28 


0.32 
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7- 


-8 
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5 


0.075 


-0.59 


± 
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3.31 


± 
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1 


.50 
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7- 
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5 
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-1.31 


± 
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2.93 


± 
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1 


.50 
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MC81 v ± 
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MHD 


7- 


-8 


3.5 


1 


0.225 


-1.00 


± 


0.07 


3.41 


± 


0.09 


1 


.18 


0.27 


MC81 t>|| 
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MHD 


7- 


-8 


3.5 


1 


0.225 


-1.09 


± 


0.06 


3.38 


± 


0.08 


1 


.18 


0.28 


ME21 v ± 
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MHD 


1- 


2 


14.0 


1 


0.5 


-0.85 


± 


0.06 


3.22 


± 


0.09 


1 


.90 


0.77 


ME21 v\\ 


128 


MHD 


1 


2 


14.0 


1 


0.5 


-0.81 


± 
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3.15 


± 


0.08 


1 


.90 


1.27 


D 
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HD/Decay 


1 


-8 


5.0 





0.1 


-0.76 


± 


0.03 


4.39 


± 


0.08 





.68 


0.19 


D 
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1 


-8 
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0.15 


-0.59 


± 


0.02 


4.70 


± 


0.08 





.54 


0.17 


D 
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1 


-8 


5.0 





0.3 


-0.32 


± 


0.02 


5.13 


± 


0.09 





.37 


0.13 


D 
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1 


-8 


5.0 





0.5 


-0.14 


± 


0.01 


5.44 


± 


0.09 





.26 


0.09 


U 
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HD/Decay 


1 


-8 


50.0 





0.1 


-0.97 


± 


0.03 


4.09 


± 


0.07 





.97 


0.35 


U 
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1 


-8 


50.0 
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-0.86 


± 


0.02 


4.48 


± 


0.09 





.81 


0.30 


u 
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1 


-8 
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-0.48 


± 


0.02 


4.93 


± 
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.51 


0.21 


u 


256 


1 


8 


50.0 





0.5 


-0.29 


± 


0.01 


5.17 




0.08 





.40 


0.15 



a Numbcr of pixels in each dimension 
b Driving wavenumber 

c rms Mach number (equilibrium M for driven turbulence; initial M for decaying turbulence) 
d Ratio of Alfven speed to sound speed 
c Timc step, measured in sound crossing times 
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Fig. 3. — Variation of P v wgt /P v true with fc for the HD and MHD fields. Fits (heavy lines) to derive <Sk were made in the range of k 
delimited by the dashed lines (largest driving scale to k = 32). 
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Fig. 4. — Variation of P wgt /P true with k for the decaying turbulence simulations as time progresses. Fits (heavy lines) to derive Sk were 

"ic '« 

made in the range of k delimited by the dashed lines (fc = 10 to k = 64). Each plot is labeled by the time step in units of the sound crossing 
time. (Left: model D, initial M = 5. Right: model U, initial M = 50.) 
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t = time since driving is turned off 
(sound crossing times) 

Fig. 5. — Plot of 8k versus t = time after driving is turned off in units of the sound crossing time. For the driven runs, we have taken 
t = 0. Filled circles are driven HD; open triangles are driven MHD; open squares are decaying HD (model D, initial M = 5); open circles are 
decaying HD (model U, initial M = 50). 
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Fig. 6. — a) Plot of 8k derived from fits to the power spectrum ratio fitted up to k = 12 versus 5ft derived from fits to the power spectrum 
ratio fitted up to k = 32, for a subset of the driven runs (see text). The solid line is a line of equality, (b) Difference between the 8k values 
shown in (fit to k = 32 minus fit to k = 12). The solid line marks a difference of zero. 
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Fig. 8. — (a) Plot of the density-weighted velocity line centroid spectral index, K2D,w, versus the unweighted (true) velocity line centroid 
spectral index, K2D u f° r an driven runs, with the indices fitted to k = 32. (b) Plot of 8k (= K2D ,u>~ K 2D ,u) versus k-id u for all driven 
runs, with the indices fitted to k = 32. (c) Plot of the density-weighted velocity line centroid spectral index, K2D,w, versus the unweighted 
(true) velocity line centroid spectral index, K2D,u for a subset of the driven runs (see text), with the indices fitted to k = 12. (b) Plot of 8k 
(= K 2D,u>~ K 2D u) versus K20 u f° r a subset of the driven runs (see text), with the indices fitted to k = 12. In all plots the solid line denotes 
the relationship 8k = -0.96 and the dashed lines mark uncertainties of ± 0.22 around this relationship. 
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Fig. 9. — Plot of 5k versus ft c( j for all simulations. Filled circles are driven HD; open triangles are driven MHD; open squares are decaying 
HD (model D, initial M = 5); open circles are decaying HD (model U, initial M = 50). The arrow shows the direction of evolution of 5k as 
time (since driving is turned off) increases; see Figure 5. 
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Fig. 10. — (a) Plot of 5k versus gt;„( p ) for all simulations, (b) Plot of re c( j versus <7;„( p ) for all simulations. Filled circles are driven HD; open 
triangles are driven MHD; open squares are decaying HD (model D, initial M = 5); open circles are decaying HD (model U, initial M = 50). 
The phase-randomized density model has cr in ( p ) = 0.22, 8k = -0.18, k c£ j = 2.87. 
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Fig. 11. — (a)Plot of o"; n (jy) versus <Ji n ( p ) for all simulations, (b) Plot of cr N / Ng versus o" p / po for all simulations. Filled circles are driven 
HD; open triangles are driven MHD; open squares are decaying HD (model D, initial M = 5); open circles are decaying HD (model U, initial 
M = 50). 
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Fig. 12. — (a) Plot of Sk versus cr p / pa for all simulations, (b) Plot of k c( j versus o"p/ po for all simulations, (c) Plot of Sk versus cf n / Nq for 
all simulations, (d) Plot of K c d versus cTjv/jv for all simulations. Filled circles arc driven HD; open triangles arc driven MHD; open squares 
are decaying HD (model D, initial M = 5); open circles are decaying HD (model U, initial M = 50). The phase-randomized model has 
°p/po = 0.21, ctjv/jv = 0.065, 5k = -0.18, and K cd = 2.87. 
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Fig. 13. — Overview of the region near NGC 7129 selected for velocity line centroid analysis. Left : 12 CO peak temperature image. Right 
13 CO peak temperature image. The wedges show antenna temperature in K. 
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Fig. 15. — (a) Comparison of the velocity line centroids derived from 12 CO and 13 CO where both measurements can be made, (b) Histogram 
of the differences in measured Vi sr ( 12 CO velocity minus 13 CO velocity). 
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Fig. 16. — The observed velocity line centroid maps. 
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Fig. 17. — Overview of the relevant terms necessary to construct the output line centroid power spectra, (a) 12 CO field; raw observed 
lower spectrum, noise floor and beam pattern (45"); (b) same for the hybrid field; the dashed line is the 46" beam pattern applicable to the 
3 CO data; no uncertainty arises from applying the 45" beam pattern to the hybrid field. 




Fig. 18. — Output velocity line centroid power spectra after noise floor subtraction and beam pattern division (Equation 5). The heavy 
lines are the fits of k^d to these spectra in the range k = 2 to k = 13. 
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Fig. 19. — Principal component analysis results for the selected field. The solid lines are the fits to the 12 CO (upper) and 13 CO (lower) 5v, 
I pairs (5v oc l a ). 
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Fig. 20. — Output integrated intensity power spectra after noise floor subtraction and beam pattern division. The heavy lines are the fits 
of K cd to these spectra in the range k = 2 to k = 13. 



